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The early stages of drop impact onto a solid surface are considered. Detailed numerical 
simulations and detailed asymptotic analysis of the process reveal a self-similar structure 
both for the velocity field and the pressure field. The latter is shown to exhibit a maxi¬ 
mum not near the impact point, but rather at the contact line. The motion of the contact 
line is furthermore shown to exhibit a ’tank treading’ motion. These observations are ap¬ 
prehended at the light of a variant of Wagner theory for liquid impact. This framework 
offers a simple analogy where the fluid motion within the impacting drop may be viewed 
as the flow induced by a flat rising expanding disk. The theoretical predictions are found 
to be in very close agreement both qualitatively and quantitatively with the numeri¬ 
cal observations for about three decades in time. Interestingly the inviscid self-similar 
impact pressure and velocities are shown to depend solely on the self-similar variables 
(r/y/i,z/y/t). The structure of the boundary layer developing along the wet substrate 
is investigated as well, and is proven to be formally analogous to that of the boundary 
layer growing in the trail of a shockwave. Interestingly, the boundary layer structure only 
depends on the impact self-similar variables. This allows to construct a seamless uniform 
analytical solution encompassing both impact and viscous effects. The depiction of the 
different dynamical fields enables to quantitatively predict observables of interest, such 
as the evolution of the integral viscous shearing force and of the net normal force. 


1. Introduction 

The impact of a liquid drop onto a rigid surface results in a rapid sequence of events 


ending, in the inertial limit, in spreading (Eggers et al. 


Hadfield|1981 1, interface tearing (Villermaux fe Bossa|2011l and ultimate fragmentation 


20101 or splashing (Stow & 


(Stow fe Stainer||1977 1. A large number of studies have investigated the many facets of 


drop impact, with a special attention to the description of its late stages (Rein 1993 


Yarin & Weiss 1995). The literature on the early stages of impact is however scarce 


in comparison. Detailed experimental data depicting the instants following impact can 


nonetheless be found in the work of Rioboo et al. (2002), that evidenced a “kinematic 


phase” where the drop merely resembles a truncated sphere and spreads as the square- 
root of time. This phase precedes the apparition of the liquid lamella. 

Probably one of the first depiction of the very first instants of drop impact dates 


back to Engel (1955). With the help of high-speed cinematography, Engel captured the 


chronology of events triggered by drop impact. He noted in particular the unvarying shape 
of the drop apex during the earliest moments of impact, which might be surprising due to 
the incompressible character of the liquid. [Engel| put forward the possible roles of inertia, 
viscosity or surface tension to explain this observation. Actually, the physical mechanism 
underpinning this behaviour is best illustrated with Figure There, the numerically 
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Figure 1. Close-ups of increasing magnitude on the pressure field developing inside an impacting 
drop in the inertial limit. The pressure field is extracted from Navier-Stokes gerris computations 
of a drop impacting a solid surface at early times (note that the surrounding gas dynamics is 
computed as well, but not represented here). Noticeably the motion is essentially pressureless 
(and therefore corresponds to a free fall) except in a concentrated region in the contact zone. 
The successive close-ups on pressure field structure in the contact region reveal a pressure peak 
near the contact line (the physical parameters are here Re = 5000, We = 250, tU/R=4:X 10~ 4 . 
The total size of the numerical axisymmetric domain is 2 R x 2 R, and the adaptive mesh has 
locally a mesh density corresponding to 32768x32768 grid points). 


computed pressure field within an impacting drop is represented shortly after impact 
(details to follow in the paper). It is readily seen that the structure of the pressure field is 
extremely concentrated near the contact zone, as in Hertz’ classic elastic contact problem. 
Conversely the outer region is essentially pressureless. This strong inhomegeneity in the 
pressure distribution therefore explains why, in the absence of any pressure hindrance, 
the upper part of the drop freely falls even after impact while remaining undeformed. 
The pressure concentration in the early stages of impact was first identified by|Josserand| 


& Zaleski (2003). From the key remark that the extent of the pressure concentration zone 


scales with the contact radius, these authors conjectured a self-similar structure for the 
pressure field and evolution with time as 1 /\ft - an hypothesis comforted by numerical 
results. Though sufficient to detect hints of self-similarity, numerical simulations were 
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Figure 2. Top: Time sequence of the pressure field developing inside an impacting drop 
(Navier-Stokes Qerris computation, fixed spatial magnification). Bottom: Corresponding trace 
of the pressure exerted by the drop on the solid substrate. The physical parameters for this 
simulation are Re = 5000 and We = 250. The snapshots correspond respectively to times 
tU/R = 10 -4 ,10 -3 and 10 -2 . 


nonetheless unable to reveal the inner structure of the contact zone until recently, es¬ 
sentially because of the very large scale ratio between this zone and the drop size. The 
increase of computational performance along with the development of adaptive numerical 
techniques for two-phase flows (IPopinet 2009) now allow to unravel the intimate structure 


fPopii 
:s[l]),C j 


of the contact zone, see Figures[1J>,c and[2] These snapshots reveal a quite complex struc¬ 
ture for the pressure, which counter-intuitively exhibits sharp maxima near the contact 
line, and not on the axis as in steady stagnation point flows. Interestingly this structure 
is reminiscent of typical pressure field structures observed in the water entry of solid 
objects, and evidenced by Wagner in the context of alighting seaplanes (Wagner 1932). 
In such problems a solid object impacts a flat liquid surface at a given velocity. Drop 
impact may be viewed as water entry’s opposite, for a liquid object impacts a rigid flat 
surface at a given velocity (see Fig. [3|. It is therefore likely that the analytical tech¬ 
niques developed since the thirties to describe with great precision the flow generated 
with the impact of an object, and proven to be in close agreement with experimental data 
(Howison et al. 19911, could be transposed for the drop impact problem. And indeed, 
building up on the analogy between water entry and drop impact, Howison et al. (20051 
proposed a theoretical investigation of two-dimensional drop impact on a thin fluid layer 
and described the different regions and scalings of importance for the flow dynamics. 
In particular, they reveal the radius of contact between the two liquids as a key length 
scaling the problem, analogously to the problem of water entry where the wet length 
of the solid is also determining and consistently with the observations of Josserand &; 


Zaleski (2003). 


The central motivation of the present study is to revisit the problem of a single spher¬ 
ical drop impacting a smooth flat solid surface at early times at the light of Wagner’s 
theory of impact, understand the dynamic fields structure and elucidate the short-time 
self-similar behaviour discerned in earlier studies. To develop a consistent theory, the 
approach followed throughout the manuscript will be to confront and cross-test sys¬ 
tematically the theoretical predictions with detailed and accurate numerical simulations 
performed with a Navier-Stokes multiphase flow solver. As a side note, we essay to make 
the paper self-contained whenever possible. In §2 we formulate the hypotheses and theo¬ 
retical framework of the problem, and describe the short-time drop impact dynamics in 
the context of Wagner’s theory. We put forward in particular a so-called “Lamb analogy” 
mirroring the flow within the impacting drop with the one induced with a flat rising ex¬ 
panding disk. In §3 we demonstrate that the Wagner flow can be recast as a self-similar 
solution for the drop impact problem. The nature of the near-axis stagnation flow and of 
the near contact line flow and pressure maxima are also discussed. Numerical results ob¬ 
tained with Cjerris (Navier-Stokes solver, VOF, adaptive mesh) taking into account surface 
tension, surrounding gas and viscous effects are compared with the theoretical prediction. 
The structure of the boundary layer is examined in §4 and is found to be reminiscent 
of the viscous boundary layer leaved in the trail of a shockwave (Mirels analogy). The 
inviscicl Wagner flow and this viscous boundary layer are found to depend on the same 
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Figure 3. Sketch of the impacting drop before contact (left), and shortly after impact (middle). 
The shape the drop would assume in absence of wall is outlined with a dashed line, and the 
contact line position is marked with red dots. This problem may be viewed as the dual of the 
classic water entry of a solid object (right). 


self-similar variable. This allows us to build a uniform, seamless solution encompassing 
both impact and viscous effects. We conclude in §5 by discussing the obtained results 
and the limits of the present investigation (such as the role of air), and by reminding 
observables of interest, such as the net impacting force of the total viscous shearing force 
exerted by an impacting drop. 


2. Model 

2.1. Theoretical framework & hypotheses 

We consider throughout this study an idealized drop impact where a perfectly spherical 
liquid drop collides with a flat rigid surface. Though classic, this model situation relies on 
a number of physical hypotheses detailed in the following. Starting with the initial perfect 
spherical shape assumption, we may identify several typical causes for deviations from 


sphericity such as capillary drop oscillations during free fall (Engel 1955 Thoroddsen 


et al. 20051 or flowing air shaping (Pruppacher & Beard 19701. Such effects will be 


disregarded in the following even though they merely result in a local curvature radius 
modification in the contact region, hence could easily be incorporated in the discussion. 
The impact flow evolution will be considered incompressible. However, the question of 
the role of compressible effects within the liquid should not be eluded when considering 
impact phenomena. A large body of literature has been devoted to acoustic effects within 
(typically fast) impacting drops. Typically, such effects are considered to play a significant 
role over lengthscales of order UR/c during timescales of order UR/ c 2 , where U stands for 
the impact velocity, R for the drop radius and c for the celerity of compressive sound waves 
in the liquid ((y) 2 C y < 1). Considering a millimetric drop impacting at a velocity 
of the order of 10 m.s^ 1 representative of e.g. a raindrop falling at terminal velocity, we 
reckon that acoustic effects matter only in a micron-sized region over a few nanoseconds 


(see Weiss & Yarin 1999 for a discussion and references). The following discussion will 
therefore be limited to those cases where the impact velocity is much lower than the speed 
of sound, as the falling raindrop, where acoustic effects can harmlessly be neglected and 
an incompressible description remains accurate. The high pressure and stresses generated 
upon impact can result in marked erosion or yielding ( Rein|199 3). Furthermore substrate 
deformation has recently been shown to significantly alter drop impact in the limit of 
very soft ( Mangili et a7||2012 l or very flexible substrates ( Antkowiak et a7]|2011 ). None 
of these effects will be considered in the following, yet an estimation of the net force 
exerted by the impacting drop on the underlying substrate will be provided in ^3] Another 
phenomenon potentially responsible for the cushioning of the impact is the thin air layer 
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between the drop and the substrate just before impact. Due to lubrication effects, this 
air layer pressurizes and dimples the drop, eventually yielding a tiny entrapped gas 
bubble in the drop (Thoroddsen et al. 2005). This phenomenon along with other roles 
of surrounding gas in impact dynamics will be disregarded in the forthcoming analysis 
(except explicitly specified, the starting state for each of the simulations is a drop already 
touching the ground - see appendix for details) and discussed in the last section. Finally 
the core hypothesis of the present study is the inertia-dominated character of impact. In 
particular, we assume that gravity, capillary and viscous effects are small with respect to 
inertial ones, i.e. Froude Fr = U 2 /gR , Weber We = pU 2 R/a and Reynolds Re = pUR/p 
numbers are all large with respect to unity. Here g denotes the gravity, cr the liquid-gas 
surface tension, p the liquid density and p its viscosity. These assumptions underpin 
the choice a purely inertial description free of these effects in the following. However, 
locally these phenomena might become more important or even dominant, examples 
being viscosity near the boundaries or capillarity in high-curvature region. In section §4 
we will address viscous effects and develop a boundary-layer correction to the inviscid 
solution, and capillary effects will eventually be discussed in §5. 


2.2. Governing equations and analogy with the water entry problem 
2.2.1. Problem statement 


We consider a perfectly spherical liquid drop of radius R and density p impacting nor¬ 
mally a flat rigid ground with velocity U, see Fig. [3j Neglecting for now the development 
of viscous rotational boundary layers, we assume that the fluid motion following impact 
is irrotational, axisymmetric and can be described with the scalar potential <p(r,z), i.e. 
the fluid velocity u(r, z ) satisfies it(r, z) = Vb>(r, z). Incompressibility requires <j> to be an 
harmonic potential satisfying Laplace’s equation, here written in cylindrical coordinates: 


ld_ ( df\ d 2 ^ 

r dr \ dr ) ^ dz 2 


( 2 . 1 ) 


The liquid dynamics obeys the unsteady form of Bernoulli’s conservation equation: 

777 + ^ I V< ^| 2 + “ = const- (2-2) 

dt 2 p 

This set of equations is completed by appropriate boundary equations. At the wall z = 0, 
the condition of impermeability reads 

^ = 0 for 0 < r < d(t), (2.3) 

where d{t) stands for the contact line position, an unknown of the problem. The position 
of the free surface is tracked with the kinematic condition: 


d S 
d t 


= 0 , 


(2.4) 


where S(r,z,t) is a function vanishing on the free surface. Expressing normal stress 
continuity at this interface yields the following dynamic boundary condition: 


p = 0 at the free surface. 


(2.5) 


Note that atmospheric pressure as here been taken as the reference pressure. 
Anticipating the forthcoming analysis of the contact region, we further note that the 
free fall behaviour outside the contact region can be recast into the following far-held 
condition: 


cj> -P- —Uz far from the contact region. 


(2.6) 
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Figure 4. Flow in an impacting drop in the fixed frame and in the drop frame computed 
with Qerris. Left, top to bottom : Streamlines and pressure map within an impacting drop for 
Re = 5000 and We = 250 in the laboratory frame at different post-impact times (t — 5 x 10~ 4 , 
5 x 10~ 4 and ICR 2 ). The overall velocity field resembles a stagnation point flow in a near-wall 
region whose extent is scaling with the wet area, and a uniform downwards flow outside. Right: 
Same velocity field in a reference frame moving with the drop initial velocity, evidencing a bypass 
motion near the contact line and an overacceleration of the free surface towards the wall. 


This condition allows to identify the constant in (2.21 as ^ U z . 

Now nondimensionalising the problem using the inertial scales R , p and U, we introduce 
the following quantities: 


R 


r = Rr , z = Rz , t=—t, (f>=UR4>, p = pU 2 p, 

and rewrite the equations into their dimensionless counterparts: 


d^ 1 (j) 

r— ( ) + -j—Tj- = 0 in the liquid, 


ld_ 

r dr 

dcj> 

sF 


,d(j) 

dr 

1 

2 


|V(/>| 2 + p= - in the liquid, 


dcj) 


dz 


(r, z = 0, t) = 0 over the wet area r < d(t), 
p = 0 on the free surface, 
d5 


df 


= 0 on the free surface. 


(2.7) 

( 2 . 8 ) 

(2.9) 

( 2 . 10 ) 

( 2 . 11 ) 

( 2 . 12 ) 


Finally the nondimensional far-held condition reads: 


= —z far from the contact region. (2-13) 

As posed, the problem entirely depends on the wet area extent d(t), whose dynamics 
has still to be determined. In the following, we investigate the near-contact line flow to 
clarify this wetting dynamics. 


2.2.2. Contact line motion: numerical observations 

To shed light on the contact line dynamics, detailed numerical simulations of impact¬ 
ing drops were carried out with Qtrris (see (j6]). Figure [4] represents typical streamlines 
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Figure 5. Advancing contact line shortly after impact. In the earliest moments following impact, 
the motion of the free surface near the contact zone is essentially directed downwards. The 
sketches show the position of the contact line for two successive instants, and illustrate the 
fact that the horizontal extension of the wet area is governed by the vertical movement of the 
interface. 


extracted from the simulations, shortly after impact. On the left panel it can be seen 
that the motion within the impacting drop far from the contact zone is vertical, uniform 
and pointing downwards, corresponding merely to the free flight behaviour —Ue z . Near 
the wall though, the flow is deflected and exhibits a stagnation point-like structure, in 
a region whose extent scales with the wet area. To investigate further the nature of this 
corrective flow, we represent on the right panel of Fig. [4] the streamlines in a reference 
frame moving with the initial velocity of the drop. There it appears that the flow winds 
around the contact line, revealing that (i) the liquid near the contact line falls faster 
than free-flight and (ii) rather than being pushed by a sweeping motion, the contact line 
progresses via a tank-treading movement, analogous to the rolling motion evidenced in 
previous studies of advancing contact lines ( Dussan V. fe Davis|1974| |Chen etHl ||1997 1. 

These observations therefore suggest that the kinematics of horizontal extension for the 
wet radius is controlled by the vertical motion of the free surface. Figure|5]illustrates this 
process, and indicates that the law of motion of the contact line d(t) can be obtained from 
the knowledge of velocity field at the free surface. Such a kinematic condition expressing 
the contact between a liquid surface and a solid object has actually been used in the 
context of the water entry of solid objects for about 80 years, and is currently referred to 
as Wagner condition. In the following we depict the analogies between these two liquid 
impact problems, and use them to derive a simple fluid mechanical model for the drop 
impact. 


2.2.3. Analogy with the water entry problem 

The modern understanding of the liquid motion and forces generated by an impacting 


object in water originates in the pioneering work of Wagner in the early thirties (Wagner 


1932). The primary motivation of Wagner was to provide a detailed characterization of 


the impulsive forces generated with impact - already known to be of sufficient amplitude 
to induce bouncing (ricochet), and even possibly structural failure of alighting seaplanes 
or slammed ships (Nethercote et al. [I~98 6). The foremost issue in this problem evidently 
stems from its highly unsteady and nonlinear nature. The central idea of Wagner was to 
model the flow induced by the impact of a float or a keel by the one induced by a flat 
“plate”, propelling the fluid particles downwards at the float or keel velocity, and having 
an extent growing with time as the waterline length. The corresponding flow (“gleiche 
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Figure 6. In the reference frame of the falling drop, the flow induced by impact may be seen as 
the one induced by a flat rising disk (Lamb disk analogy). The winding motion is here represented 
with orange arrows, and the radial expansion of the disk with the wet area is indicated with 
purple arrows. The motion of the disk itself is given by the red arrow. 


Tragfiugelbewegung” - equivalent aerofoil motion) is typically found to wind around the 
plate and therefore to promote jetting or splashing. The knowledge of this flow field 
then allows to determine the motion of the free surface, and finally provides the needed 
condition in the determination of the wet length d(t). 

Analogously, for the drop impact problem, our numerical simulations evidence similar 
flow features and winding motion. These observations advocate for the use of a water 
entry-analogue description, where the flow induced by drop impact would correspond to 
the one induced by a flat expanding disk in normal incidence, which extent is given the 
wet area (see figure [b]). Following this vision of drop impact as a dual version of the water 
entry problem, we adopt from now on the corresponding formalism to describe the fluid 
mechanics of impact. 


2.3. Leading-order description for the drop impact problem 

Interested in the early-time behaviour of the impact-induced flow, we set out by exam¬ 
ining time-dependent solutions of system (2.8 - 2.12) in the vicinity of the contact zone. 
To this end, we introduce e as a measure of the wet region: d(t)/R = 0(e) (see Fig. [7]). 
This £ is the fundamental small parameter of our problem. 

As typical in two-phase phenomena, the lengthscales for the dynamical fields and for 
the geometry of the free surface differ in this problem. Starting by considering the space 
variables f and z on which depend the dynamical fields (such as the velocity potential 
</> or the pressure p), we introduce the following rescaling: r = e r f and z = e~z, where 
r and z are 0(1) quantities and e r and e z are gauge functions. From the structure of 
Laplace operator, we expect the dynamical fields to display identical length scales in each 
direction, so that e r = £ z = £■ 

Insights into the relevant lengthscales for the description of the free surface geometry 
can be gained by decomposing the position of the surface into that of a translating 
sphere Zs{f,t ) plus a surface disturbance h(r,t ) (see Fig. !>)■ Assuming the drop falls 
with constant velocity, the shape of the unperturbed translating sphere obeys r 2 + (zs — 
(1 — fj) 2 = 1. Sufficiently close to the contact area, we introduce gauge functions for 
the vertical position of the moving sphere zs and the time t: zs = £ zs zs and t = syt. 
The equation for the sphere surface can be approximated by e zs zs = |e 2 r 2 ~ e tt- As 
previously the determination of these scaling functions is obtained by dominant balance 
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♦ £ zs = e 


£ r = £ 


£ r = £ 


Figure 7. Scalings in the contact zone. At the earliest times only a very small portion (of order 
e) of the drop touches the wall. The fluid sets into motion with impact is in a region of extent e 
in every direction. The air wedge confined between the wall and the drop presents an angle of 
order e as well. The colormap illustrates the pressure distribution. The physical parameters for 
this simulation are Re = 5000 and We = 250. This snapshot corresponds to a nondimensional 
time f = 10 -3 . The position of the contact line is here d= 1.73 10~ 3 ^ 2 . 


arguments: e zs = = e 2 . Note that at short times the intersection radius between the 

sphere and the impacting plane is given by fj n t er sect = 

We remark that as in the original study of Wagner, a scale separation between zs and z 
exists ( small deadrise angle hypothesis , see e.g. Qliver|2002 ). This scale separation arises 
because the drop typical radius of curvature (0(1)) is very large in front of the other 
lengthscales of the problem, see Fig. [7] 

We now turn on to the surface perturbation h(r,t), that embodies the impact-induced 
flow. Recalling that h represents a perturbation around a falling sphere, we can express 
the position of the free surface with the following implicit equation: S[f,z,t) = z — 
[zs — h(r, £)) . Introducing an appropriate gauge function £h such that h = e^h we obtain 
by dominant balance analysis that = e 2 . It follows that: 


S[r,z,t) = z - - f 2 + t + h[r,t) 

= 0 on the free surface. 


(2.14) 


The kinematic boundary condition derives from the previous equation. At the free surface, 
we have: 

dh d(j> 


AS dh _ df) 

dF = 

0 [ 1) O(e 0 ) 


■ £cf> 


dr dr 

ol ^7 


_ n 

' JJl -°> 

0{eje) 


(2.15) 


where 0 = e^0. It is impossible here to keep all terms at the same order; the dominant 
balance between the vertical velocities dh/dt and d(j)/dz implies that £$ = e. At leading 
order, the kinematic boundary condition is therefore reduced to : 


(90 

dz 


dh 

~dt 


1+-—; + —^ = 0 on the free surface. 


(2.16) 


It proves convenient to introduce a translation of the velocity potential such that 0 = 
—z + 0. This translation merely accounts to analyse the problem in the falling-drop 
reference frame. The kinematic boundary condition is then simply rewritten as: 


dh 

~dt 


(90 

dz 


on the free surface. 


(2.17) 
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r = er z = ez t = e 2 t 

p = s~ 1 p (j = ecj> (u, v) = (u, v). 

Table 1. Resume of the most important asymptotic scales of the problem. 


Inserting these different scaled variables into Bernoulli’s equation, we obtain: 


e p p- 


1 d<j) 

e dt 



= - in the liquid, 


(2.18) 


where p = £ p p. The scale of the pressure, £ p = —, is here seen to be as large as the contact 
zone is small - as expected in an impact problem. At leading order, Bernoulli’s equation 
is therefore reduced to: 

j n j.jjg liquid (2.19) 

at 

It follows from this equation that the constant pressure Dirichlet boundary condition 
on the free surface p = 0 can be recast as a condition for the potential at the free surface: 
(f> = const, where the constant is arbitrary. Without loss of generality, we set from now 
on this constant to zero. 

Finally, as classic in water wave theory, we exploit the shallowness of the gap between 
the free surface and the plane to transfer the boundary condition at the free surface onto 
the plane (see e.g. Van Dyke||1975 §3.8). 

Summarizing, the model problem takes the following expression: 


1 d ( dcj)\ d 2 (j) 

--^z r— + wvy = 0 m the liquid, 
r or y or I oz z 

d(j) 

— —- = p m the liquid, 
at 

the locus d(t) of the contact line is determined with the Wagner condition: 

h{r,t) = \r 2 -> 
so that the boundary conditions at z = 0 read: 


( 2 . 20 ) 

( 2 . 21 ) 


dh 


and the far-fielcl behaviour is given by: 


dtp 
dz 

= , 
dz 


II 

(2.22) 

for r > d(t), 

(2.23) 

for f > d(t), 

(2.24) 

for r < d(t ), 

(2.25) 

as r, z —> oo 

(2.26) 

as f —> oo. 

(2.27) 


Finally the corresponding model geometry is sketched Fig. [8] We remark that the previous 
set of equations resembles to that of the classic water entry problem, and can be solved 
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A(j> = 0 


<t> = 0 <t> = 0 

-i-1- 

—d(t) +d(t ) 

Figure 8. Leading order outer problem for times of order e 2 . 


using the methodology described in e.g. Oliver (2002). In the next section though we will 
present an alternate method based on self-similar solutions. 


3. Self-similar solutions and numerical simulations 

3.1. A self-similar problem 

To reveal the self-similar nature of our problem, we classically look in the following for 


scale invariance (Darrozes & Frangois (1982)). We start by expressing the fact that any 
variable q in (f, z, t, f>, h, d,p) can be rewritten as q = X q q, where q is a rescaled variable 
and X q a numerical stretching coefficient embodying the change of scale. Inserting these 
variables into the governing equations, it is straightforward to see that invariance of 
Laplace equation through this stretching requires X r = X z . Similarly, expressing the 
invariance of Wagner condition yields Xh = Xt, X r = y/Xt and A d = y/Xt- The same 
operation performed on the additional boundary conditions finally imposes A^ = VV 
and X p = 1/y/Xt. Note that A* remains here as the sole stretching parameter. 

The pressure field can be written as an implicit function of time and space as fol¬ 
lows: F{p,f,z,t) = 0. Upon using the previous scale invariance arguments, this rela¬ 
tion may be rewritten as F{p/y/Xt, y/X/r, \fXt.z, Xtt) = 0 . A simple algebraic manip¬ 
ulation allows to remove the A t dependence for all but one variables, so that finally 
G(Vtp, r /Vi, z/Vi, X t t) = 0 , for any A t . Remarking that for a given t, this function has 
to cancel whatever the choice of the scale A*, it readily appears that the last variable is 
superfluous. In other words, a relation linking Vt p to f/Vt and z/Vt only must exist. 
The pressure field may therefore be rewritten explicitly as: 


1 / f z 

p= vf\7tTt 


(3.1) 


With a similar reasoning, and upon introducing the self-similar variables f = f/Vt and 
V = z/Vt, we readily obtain : 

(f(r,z,t) = \/i <!>(£, 77 ), h(r,t) = t TL(f,) and d{x,i) = yft S, (3.2) 


where *I> and T~L are unknown functions of the self-similar variables and S a constant 
representing the (fixed) position of the contact line in self-similar space. This allows us 
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to formulate the self-similar version of the drop impact problem : 


19 / 9<f> 

£d£ 


<9 2 <f> 

dij 2 


= 0 in the liquid, 


9$ 9$ 


V{£,r}) = - ^-$(£,77) + +rt ~d^J in the lic l uid ’ 
the boundary conditions at 77 = 0 take the following form: 

1 on 9 $ 

n -~2 for {>i - 
9 $ 

~drj = 1 fOT *<*’ 

$ = 0 for £ > 6 , 

the far-field behaviour is: 

$ —» 0 as £, 77 —> 00 
H —> 0 as £ — > 00, 

and the self-similar version of Wagner condition is finally given by: 

^(0 = \Z 2 - 1 for £, = 5. 

This problem can now be solved in several steps. 


3.2. Self-similar potential 

In this geometry, Laplace equation can be solved with variable separation, leading to a 
family of elementary cylindrical harmonic solutions with an exponential behaviour in 77 
and an oscillatory one in £. We recompose by summation and obtain : 


(3.3) 

(3.4) 

(3.5) 

(3.6) 

(3.7) 


(3.8) 

(3.9) 


(3.10) 


poo 

*(£, 77 )=/ C(k)Jo(k£)e~ kri dk. 
Jo 


(3.11) 


The weight function C(k) is determined with boundary conditions (3.6I and (3.7), leading 
to the following pair of dual integral equations: 


pOO 

/ kC(k)J 0 (k£) dk = -1 for £ < 5, 

Jo 


C(k)J 0 (k £) dk = 0 


for £ > S. 


(3.12a) 

(3.126) 


Solving these dual integral equations using the technique described in Sneddon (1960), 
we obtain a closed-form expression for the weight function: 


2 Skcos(kS) — sin(kS) 2 d /sin(kS) 
7 r k 2 n dk \ k 


(3.13) 


Anticipating the description of the contact line dynamics, we now derive dd>/drj at the 
substrate level 77 = 0 : 


94> 2 r°° kScos(kS) — sin(/c<5) 

9t? 7T J 0 k 2 


Jo(k£)kdk, 


(3.14) 
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where we recognize the sum of two Hankel transforms (see e.g. Sneddoir| |1995[ table IV, 
page 528). This allows us to obtain the following explicit expression for /dr] for 77 = 0: 


<9$ 9$ 2 1 S 

— = 1 for £ < 0 and — =- , 

dr] dg t: l y^ 2 - S 2 


— arcsm [ - 



for £ > S. (3.15) 


3.3. Wagner condition and contact line dynamics 
With the help of the vertical velocity expression in the near wall region just derived, we 


can rewrite the kinematic boundary condition (3.51 as: 




— arcsm | - 



for £ > 5. 


(3.16) 


This inhomogeneous differential equation can be solved using variation of parameters, 
i.e. looking for a solution of the form 'H(f) = £ 2 /(£)- This gives: 


m 


+00 


J s 




(3.17) 


Upon using the far-field decaying behaviour of H (see equation (3.9)), this last equation 


reduces to f(S) = |<5 2 so that at the contact line the drop deformation is: 

H(5) = 5 2 f(5)= 1 -. 


(3.18) 


In the self-similar space, the Wagner condition therefore takes the following remarkably 
simple form: 


- = -< 5 2 — 1 , 

2 2 

from which we finally derive the position of the contact line: 

5 = V3. 


(3.19) 

(3.20) 

It is interesting to remark that the contact line motion d(t) = V3t just predicted 
within the framework of Wagner theo ry is quite close from the rough truncated sphere 
approximation fintersect = v^2 1 (Rioboo et al. [20021. Furthermore, and as for other liquid 
impact problems, we note that the agreement between this prediction and observations 
is fairly satisfactory. Indeed Fig. [9] reports early post-impact successive positions of the 
contact line extracted from numerical simulations performed with Qerris along with our 
theoretical prediction. Noticeably the superposition between theory and numerical results 
is excellent, at least until the moment of formation of a liquid corolla (here indicated with 
a red dot). 

3.4. Analogy with the normal motion of an expanding disk in an infinite mass of liquid 
In §2.2.3| we proposed to visualize the flow in an impacting drop as the one induced by a 
flat rising disk expanding radially as the wet area (see also Fig. | 6 |). We are now in a posi¬ 
tion to formally justify this water-entry analogy. The axisymmetric flow induced by ‘the 
motion of a thin circular disk with velocity U normal to its plane, in a infinite mass of liq¬ 
uid’ is for example analysed in Lamb’s classic textbook §101 (Lamb il932| . After deriving 
some elementary axisymmetric solutions of Laplace equation of the form exp(±fc^) Jo(fcr) 
in §100, [Lamb examined a variety of axisymmetric potential flows. Among those was the 
one (later connected to the flow around a flat circular disk in normal incidence) where 
at the symmetry plane z = 0 the potential takes the value </> = C\]a? — r 2 for r < a and 
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Figure 9. Comparison between the theoretical position of the contact line as a function of 
time deduced from Wagner theory d(t) = v3 1 (dashed line) and the position of the contact line 
extracted from Cjerris computations of an impacting drop at Re = 5000 and We = 250 (grey 
dots). The red dot marks the birth of the corolla. 


(f> = 0 for r > a, with a the disk radius. The solution for this problem was stated under 
the following integral representation: 

<Kr, z) = -C j™ e~ kz Mkr ) A d k. (3.21) 

And from ‘a known theorem in Electrostatics’, Lamb obtained the expression for the 
vertical velocity in the symmetry plane: 

-7 tC for r < a, (3.22a) 

C ^arcsin ^ ^ for r > a. (3.226) 

This corresponds precisely to the flow within the impacting drop, after posing C = —2/7r 
and a = S, thereby justifying formally our initial analogy between the impact-induced 
flow with the one associated with a flat rising disk rapidly expanding with the wet area. 
Setting C = ZU/tt, Lamb remarked that the above potential indeed describes the flow 
winding around a flat disk moving at velocity U. He further noted that a simple expression 
for the fluid half-space kinetic energy could be derived from the previous relation: 

Tdisk = ^pa 3 U 2 . (3.23) 

This expression can immediately be transposed into the (nondimensional) kinetic energy 
of the impact-induced flow within the drop: 

T = AV3P /2 , (3.24) 

or, equivalently, into its dimensioned counterpart: 

T = 4V3pU 7/2 R 3/2 t 3/2 . (3.25) 

We emphasize that this expression is derived within the frame of the falling drop and, as 
such, represents the kinetic energy of the defect flow associated with impact. Although a 
direct physical interpretation of this quantity is not straightforward, we will see in §3.7| 
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that the knowledge of this defect kinetic energy will allow for a direct determination of 
the impacting force. 


3.5. Structure of the velocity field 

We now investigate the structure of the velocity field in the contact region and search 
for exact closed-form expressions and convenient approximations for this field. 


3.5.1. Integral representation of the velocity field 

In the fixed frame, the velocity field u(f, 5, t) inside the impacting drop can formally be 
derived from the (untranslated) potential —77 + d>. Following the arguments developed in 
S |3.1| this velocity field is simply related to the self-similar velocity field U{^rf) via the 
relation: 

u(r,z,t)=U(£,rj). (3.26) 

where the components of the self-similar vector field IA = {U^. U rl ) are: 


! <9<F 


(3.27a) 

(3.276) 


Inserting the expression of the self-similar potential determined previously yields the 
following integral representation for the vector field components: 


2 r°° 

U^,rf) = -/ 

n Jo 

W„(£,?7) = -1 - - 

7r 


■Akc^k )-si»('/^) e - blJi(fa;) 
Jc 

In h 


(3.28 a) 
(3.286) 


A closed-form expression is unfortunately not accessible in the general case. In the fol¬ 
lowing however we calculate the value of these integrals at some particular places. 


3.5.2. Closed-formed expressions for the velocity field along the axis and the substrate 


Simple analytical solutions for the velocity field can be obtained from (3.28) at precise 


locations. Along the symmetry axis for example, where £ = 0, the properties of integrals 
of exponentials allow to write: 


W„(£ = 0,7?) = -l + 




arctan 


V3 

V 


V3y 

3 + r] 2 


for 77 > 0 . 


(3.29) 


This last result is confronted with numerical velocity profiles extracted from Qzrris compu¬ 
tations in Fig.|10| The nice agreement between the theoretical solution and the numerical 


profiles seen in the self-similar space (Fig. 10 f) here holds over more than a decade in 
time. 

Analogously, analytical forms for ( 3.28| ) can also be obtained along the substrate plane 
77 = 0 by exploiting the properties of Hankel transforms (Sneddon| |1995 1. An expression 
for the vertical velocity is already provided with equation (3.151, after inserting 5 = y/3. 


We remark that this analytical solution elucidates the faster-than-free-flight motion of 
the free surface near the contact line discerned in §2.2.2| Likewise the radial velocity 
distribution across the wet area is found to be: 


(£, 77 = 0 ) =- s 

s ^v/3^2 


for 0 < £ < V3. 


(3.30) 
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-u z (r = 0 ,z,t) 



0 0.2 0.4 0.6 0.8 1.0 

-U v (£ = 0,ri) 


Figure 10. Left : Axial velocity profiles along the axis extracted from gerris computations at 
times t = 3 x 10 -3 , 5 x 10~ 3 , 10~ J , 5 x 1CF 2 and 1CP 1 . R ight : Comparison between the analyt¬ 
ical prediction for the axial velocity given by equation (|3.29|) (red dashed line) and numerical 
solutions obtained with gerris, rescaled in the self-similar space (blue solid lines). The physical 
parameters for this simulation are Re = 5000 and We = 250. 


This unphysical inviscid slip velocity u e (tf) cannot be observed in our simulations encom¬ 
passing viscous effects. But this quantity is nonetheless relevant for it corresponds to the 
edge velocity of the viscous boundary layer (studied in detail in (jdj) . 


3.5.3. An unusual stagnation point flow 

In the very vicinity of the origin, the first order power series of the velocity field (3.28 1 
reads: 


or, equivalently, in dimensioned variables: 


u r (r, z, t) 
u z (r,z,t) 


2 fU r 

7t\/3 V R \ft ’ 
4 fU z 

7TA/3 V R \/t 


(3.31a) 

(3.316) 


(3.32a) 

(3.326) 


Though simple, this peculiar structure for the impact-induced unsteady stagnation point 
flow is nonetheless counter-intuitive and could not have been inferred from simple dimen¬ 
sional analysis. Noteworthy enough, this result is at variance with the typical structure 


of the later intermediate flow associated with spreading u ~ (r/t, —2z/t) (see e.g. Eggers 
et aLpOlO Lagubeau et al. |[20 12 Yarin fe Weiss|[l995 1. 


3.5.4. Beyond the stagnation point: a remark on the overall velocity field structure 
The previous approximation for the impact flow is valid in a small region near the origin. 
To further investigate the limits of this representation we show Fig. different radial 
velocity profiles corresponding to various locations f. The collapse of the numerical pro¬ 
files taken at various f and t (but such that r/Vt = r/Vt is constant in each figure) 
onto the theoretical profiles is again an illustration of the relevance of the self-similar 
representation. But it is also to be noted that while the stagnation point ansatz disre¬ 
gards any radial velocity variation in the profiles exhibit a sensible variation along the 
vertical coordinate 77 . This variation is best depicted with Fig. [12] where theoretical radial 
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0 
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0.4 0.6 


0.8 


(£) {=1.5 


1.0 


Figure 11. Self-similar radial velocity as a function of ri for £ = 0.125, 0.25, 0.5, 1 and 1.5. Thes e 
velocities are rescaled by the outer solution of the boundary layer w e (£) given by equation (3.30|). 
Blue solid lines represent the numerical solutions extracted from gerris computations in the 
self-similar space for t = 5 x 10 —3 ,10~ 2 , 5 x 10~ 2 and 10 _1 (Re = 5000 and We = 250). The red 
dashed line represents the theoretical solution U^[f, rj) given by equation (|3.28a|). Note that the 
boundary layer is so thin that it is almost indistinguishable (see also Fig.l 


velocity profiles taken at different values £ have been represented. Noteworthy enough, 
profiles corresponding to £ < 1 collapse on a single curve. We therefore speculate the 
following functional dependence for the velocity = u e {f ) f iji) in this region. For 

larger values of £ though, significant deviations from this behaviour arise and variable 
separation cease to hold: Uc(f,rj) = u e {09{^ r i) w ith = 0) = 1 for £ > 1. 
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u e { 0 

Figure 12. Evolution of the analytical self-similar radial velocity given by equation (|3.28a|) as 
a function of p for £ = 0.125, 0.25, 0.5, 1 and 1.5. These vel ocitie s are rescaled by the outer 
solution of the boundary layer U^{^,p = 0) = u e (£), equation (|3.30[l. 



Figure 13. Comparison between the flow pattern within an impacting drop (left) and around 
a rapidly expanding disk (Lamb analogy, right) in the self-similar space. In both cases, the 
streamlines are represented in the moving frame. The red dots represent the theoretical position 
of the contact line £ = \/3. The numerical streamlines represented on the left are derived from 
the velocity field computed with gerris at t = 10 -3 (Re = 5000 and We = 250). The theo retical 
streamlines shown on the right correspond to isovalues of '3/ (^, p) defined in equation ( 3.34|| (note 
the correspondence with Lamb’s figure page 145). 


3.5.5. Flow pattern, contact line bypass and Lamb analogy 

We now define the self-similar stream function 4) (£, 77 ) in the drop reference frame from 
the potential: 


1 d'Sf 54> 

£ drj 5£ 1 

154 ' 54 ) 

£ 5£ dp 


(3.33a) 

(3.33&) 


By integration of the previous relations we deduce the following expression for 4 f (£, 7 y): 


= - 
7r 


V3tcos(V3t)-sin(V3q e _ t , {Ji(t{) ^ (3 34) 

nj 


up to a constant. Formally, 4'(£> 7 7) is the stream function describing the winding flow 
around a flat rising disk (Lamb 1932 §108). Figure 13 offers a comparison between 


the streamlines of this Lamb analogy and the ones computed with Qaris for the drop 
impact problem in the self-similar space. A good qualitative agreement between the 
analytical and the numerical streamlines is noticeable, comforting the expanding disk 
analogy followed here. Interestingly the winding motion around the contact line, as well 
as the falling velocity overshoot near this region, are both captured with this analogy 
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Figure 14. Comparison between the pressure field developing inside an impacting drop (left) 
and around a rapidly expanding disk (Lamb analogy, right). The pressure field represented on 
the left is extracted from Qerris computations and represented in the self-similar space (t = 10“ 3 , 
Re = 5000, We = 250; The isovalues are: 0.12, 0.24, 0.36, 0.48, 0.6, 0.72, 0.84) . The self-similar 
theoretical pressure field represented on the right is given by equation ( |3.4| ) (isovalues: 0.13, 
0.28, 0.445, 0.57, 0.73, 0.9, 1.2). Though isovalues have been slightly changed between the two 
panels, theoretical and numerical results are in a good overall agreement. 

and can be correlated with the peculiarities of the winding flow near the edge of a rising 
disk. 


3.6. Self-similar pressure 

From the knowledge of the velocity potential we are now in a position to derive the 
pressure field as the time derivative of the potential. In the self-similar space, the pressure 
field is given by equation (3.4 (. Figure 14 proposes a comparison between the structure 
of the self-similar pressure extracted from numerical computations performed with Cjerris 
and the theoretical prediction. There it can be seen that the overall structure of the 
pressure field developing in the impacting drop, and in particular the pressure peak in 
the vicinity of the contact line already pinpointed out in Fig. [T] nicely matches with 
the theory. Interestingly, the structure just described is at variance with the pressure 
distribution around a flat disk rising steadily (Lamb’s original problem). Indeed in such 
a configuration the pressure is expected to be maximal in the stagnation point area, 
whereas in our model problem the pressure peaks near the contact line/disk edge. This is 
a consequence of the motion unsteadiness: the pressure is here dominated by the d(f>/dt 
contribution rather than the steady |V<^ 2 term. 

As in §3.5.2| closed form expressions for the pressure can be obtained along the axis 

and the substrate plane. The radial structure of the self-similar pressure across the wet 

area reads V(f,r] = 0) = — 2 — r for 0 < £ < x/3. This analytical prediction is con- 

tV 3 -? 2 


fronted with Qerris numerical results in Fig. 15 After a transient numerical initialization 
phase (corresponding to the red curves), the pressure profiles collapse on the self-similar 
analytical solutions (blue curves). In accordance with the overall pressure field structure 
depicted earlier, the pressure radial profile presents a local minimum at £ = 0 and a 
maximum in the vicinity of the contact line, that is for £ = \/3 - where the analyti¬ 
cal solution exhibits an inverse square-root singularity. We note that for later times the 
pressure peak is smoothed out in the numerical simulations (grey curves). As this regu¬ 
larization coincides with the birth of the ejecta sheet, we conjecture that this fall-off can 


appropriately be described with a second-order Wagner theory (Korobkin 2007 Oliver 


2007). 






















Figure 15. Left: Pressure trace on the substrate 2 = 0 obtained from the numerical simulations 
between t = 5 x 10~ 4 and 10” 1 for Re = 5000 and We = 250. The color code for each decade is 
the same as in Fig. [IT] Note that the curves are equally distributed within each decade. Right: 
Same as in the left but in the self-similar space. The black dashed curve represent the analytical 
solution for V(£,ri = 0). 



Figure 16. Pressure alon g the axis r = 0 obtained with Qe,rris for t = 5 x 10 


1 10 


" 2 ,5 x 10~ 2 


and 


10 (same as in Fig. 10 right) represented in the self-similar space (Re = 5000 and We = 250). 


The red dashed line is the analytical solution for V(£ = 0, rj). Fluctuations of the pressure around 
the theoretical prediction is to be related with numerical projections errors. 


Similarly the expression for the self-similar pressure along the s ymm etry axis can also 

3V3 


be obtained analytically: V(£ = 0 , 77 ) = 


^( 3 +^) for V > 0. Figure 16 
result with rescaled axial pressure profiles extracted from numerical" 


compares this last 
simulations. There 


again the agreement between the computations and the theory is seen to hold for a large 
time span. 

The structure of the pressure field in the vicinity of the origin can be inferred from 
these last results. Reexpressing the pressure cuts determined in terms of f, z and t, we 
get: 


f p(r,z = 0 ,t) 
[ p(r = 0 ,5, i) 


3 

7t\/3 t — r 2 
3 ^ 3 1 
7 t(3 t + z 2 ) 


(3.35a) 

(3.356) 


From these two relations, we deduce the following expansion for the pressure near the 
origin: 

«f,2,t)=^H(l + T-£_^ + .- (3 . 36 ) 

This expression provides with a local approximation for d<j>/dt from which, after time in¬ 
tegration and space differentiation, we readily recover the stagnation point flow structure 
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Figure 17. Time evolution of the pressure p(0,0,£) measured at the origin in the numerical 
simulations with gerris (Re = 5000 and We = 250). Note that each decade is represented with 
a different colour. The theoretical prediction p(0, 0,t) = ^ t _ 5 is superimposed with a black 
dashed line. 


found earlier: ( Uf,Uz ) = —2 z/Vt). This near-axis behaviour emphasizes again 

that simple intuitive dimensional analysis suggestion r/t and —z/t is here not relevant. 
The leading order term for the pressure at the origin follows: 


p(0,0,t) = — t J 
7r 


or, with dimensions p(0,0,t) = 


pU 3 ' 2 



(3.37) 


Josserand & Zaleski 


(2003) on the 


This result extends the t~ s scaling law proposed by 
basis of scaling arguments. A comparison between this theoretical prediction and Cjerris 
numerical simulations is proposed Fig. |17| using the color code of Fig. [15] After a numeri¬ 
cal transient phase, the pressure rapidly reaches the self-similar regime. Remarkably, this 
short-time similarity regime holds for almost 4 decades in time and is nicely captured 
with Wagner theory. Eventually a departure from similarity occurs when t becomes of 
order 1, and the pressure promptly drops to 0. 

3.7. Normal force induced by drop impact 

Building on the last set of results, we deduce the total net normal force imparted by an 
impacting drop on the underlying substrate at early times. Integrating the pressure on 
the wet surface, we have: 

Fit) = JJ ^ V& p = 0)dS = 2 nV~t J 3 —d£ = 6 Vat. (3.38) 

The dimensional counterpart of this net total force induced by the drop on the substrate 
therefore reads: 

F(t) = 6VSpU 5/2 R 3/2 Vt, (3.39) 

where the force is seen to increase as t? for short times. Interestingly F(t) could have been 
inferred directly from energy arguments, with no knowledge of the pressure distribution. 
Indeed, writing the global kinetic energy conservation for the upper semi-infinite space, 
we have: 

pu ■ ndS , where T = ( /// ^r- dE^ . (3.40) 


^T= - 
dt 


In the context of a flat rising disk, the kinetic energy reduces to Xdi S k = |pa 3 f7 2 (Lamb 
§102). This expression can immediately be t ransp ose d to the impacting drop problem so 
that T = dV3pU 7 / 2 R 3 Dt 3 / 2 (see equation (3.251 in §T4j). The power of pressure forces 


1 /2 

then follows as ^T = 6\/3pU 3 R 2 () . Dividing this power by U, we recover exactly 
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the previously obtained result for the net normal total force. This alternate derivation 
of the normal force provides with yet an other illustration of the relevance of Lamb’s 
analogy for the drop impact problem. 


4. Matching with the viscous solution 

The inertial limit (large Reynolds number hypothesis) investigated so far has allowed 
us to model the flow within an impacting drop as the winding motion of an inviscid fluid 
around an expanding disk, appropriately described by an harmonic potential obeying 
the unsteady Bernoulli equation 2.21. Actually the agreement between the correspond¬ 
ing theoretical results and numerical Navier-Stokes computations carried out with Cferris 
(encompassing viscous effects) comforted this approximation, see e.g. Figs 10 for veloc¬ 
ity, [16] for pressure or [9] for contact line motion comparisons. Most presumably, viscous 
effects are here dominating only in very thin boundary layers developing along the wet 
substrate. And indeed, even if the overall agreement between the radial velocity profiles 
and the inviscid solution is evident, a careful examination of Figure [Tl] reveals the pres¬ 
ence of these thin layers in the very vicinity of the solid wall. Even if spatially confined, 
these boundary layers nonetheless play a key role when comes e.g. the question of the 
erosion potential of an impacting drop. Consequently we now set out to describe the in¬ 
ner structure of these viscous layers and to match it to the previously determined outer 
inviscid solution. Viscous shear stresses and total erosion potential are eventually briefly 
discussed. 


4.1. A simple boundary layer problem? 

Typically, the (inviscid) slip velocity u e (r,t ) = ^f/y/St — r 2 , here first introduced equa¬ 
tion (3.301, and the no-slip condition at the substrate, trademark of real fluids, are 
reconciled through the introduction of a viscous boundary layer. According to the classic 
boundary layer theory (e.g. Schlichting|l968), the transverse scale of this layer is Re _1,/2 , 
so that an appropriate inner coordinate Z can be defined via z = Re - ' Z. The most 
simple idea at this point is to think that the outer variables scales defined in j )2.3| imply 
that the non-linear terms of the boundary layer equation are negligible when compared 
to unsteady and viscous terms, so that this equation would simply read: 


dU r dp d 2 U r 

dt dr dZ 2 ’ 


(4.1) 


where capitalized variables refer to boundary layer quantities. Considering that Euler 
equation in the inviscid outer domain reduces to = —1|, the boundary layer equation 
can be recast as the following diffusion equation for the defect velocity: 



u e ) 



U e ). 


(4.2) 


The corresponding solution can then be expressed as a convolution between the forcing 
term and the Green function of the heat equation: 


U r = u e (r,t) 




Ue(r,r) 


dr. 


(4.3) 


Unfortunately this solution leads to a paradoxical cancelling of shear stresses at the wall. 
We conjecture that this unreasonable result stems from the fact that the sharp longi¬ 
tudinal variations associated with the contact line motion have here been disregarded. 
Specifically non linear terms do balance unsteady terms, at least near the contact line 
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location r = V3t. As a result, the boundary layer actually grows from this moving point 
both in space and time. While a comprehensive analysis of this problem demands a care¬ 
ful balance of each term likely resulting in a non linear boundary layer problem, beyond 
the scope of the present study, we nonetheless propose in the following an approximation 
based on an analogy with boundary layers developing behind shockwaves. 


4.2. Approximation of the drop impact boundary layer via an analogy with 
shock-induced boundary layers 

We now depict qualitatively the inner viscous structure of the velocity field by using a 
simple analogy. First remembering the tank-treading movement in the vicinity of the 
contact line observed and discussed in §2.2.2[ we point out the violent change in radial 
velocity when passing through the contact line. In other words, the contact line embodies 
a neat discontinuity where the slip velocity sees its value suddenly change from 0 to u e . 
Building on this observation, we consider in the following the contact line as a kind of 
shock wave sweeping the substrate, and seeding a boundary layer in its trail (see figure 
18). This problem is classic in compressible flows and was solved by Mirels (19551 in 


the context of a shock tube (see Schlichting 1968 for more details). In this study, a 
fluid initially at rest is swept by a shockwave travelling at celerity U s in the direction 
x and instantly acquires an impulse of velocity Uoo in the process. Behind the normal 
shockwave is left a growing viscous boundary layer. 

The Ansatz for Mirel’s solution is to introduce r) m = zj\Jt — x/U s as the self-similar 
variable. This variable not only takes into account time variations but also longitudinal 
effects from the shock backwards in x. Disregarding any pressure gradient but considering 
both unsteady and nonlinear effects, the momentum equation may be rewritten in terms 
of and of the velocity Uoaf'^m): 

f'"(Vm)+^(Vm- l ^f{Vm))f ,, (Vm) = 0, with /(0) = /'(0) = 0, and /'(oo) = 1. (4.4) 


Note that compressible effects have here been absorbed via an appropriate Lees-Dorod- 


nitsyn’s transformation (see Stewartson 1964). Two limiting cases clearly emerge from 
the picture. For large Uoo/U s (and after a rescaling and a change of sign due to the choice 
of origin), the velocity profile tends to a Blasius profile. Conversely, for small values of 
the velocity ratio, the velocity rather adopts an error function profile. Note that profiles 
corresponding to intermediate values of this ratio can be found in Schlichting’s textbook. 

From this sound result we may by analogy transpose this approach to the drop impact 


problem (see Fig. 18). Obviously the outer solution for the drop impact problem is quite 
more complex for neither U s nor U Q0 are constant. The core idea consists in drawing a 
parallel between the shock (at position U s t) and the contact line (at position V / 3t) on 
the one hand, and between the steady slip velocity Uoo and u e (f, t) on the other hand. 
Following this simple analogy the longitudinal velocity is approximated with: 


U r (f,z,t) = 


2 f 


f 


- 


VRe ] . 


V3t - f 2 \2^/i-r 2 /3 ) 


(4.5) 


where f is solution of an equation which is analogous to Eq. (4.41. The so-called com¬ 
posite solution (Van Dyke[1975), which is an expansion valid in the ideal fluid and in the 
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Figure 18. Left: Sketch of the contact line during its motion and of the growing boundary 
layer in its trail, analogous to that developing behind a shockwave. Right: S hockwave-induced 
boundary layer, reproduced from the german edition of Schlichting textbook ( Schlichting| 1968 1. 
Notations are from Schlichting, with a correspondence between x and r. Note that in the shock- 
wave case, Uoa and U s are both constant. 


boundary layer, then follows: 
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2 \Jt — f 2 /3 J 


In practice we approximated f with erf function. Figure 19 proposes a comparison be¬ 
tween the numerical velocity profiles extracted from Qzrris computations and this approx¬ 
imation, which proves to provide a fairly good description for the flow. As a side note, 
we remark that replacing the error function with Blasius profile yields a slightly more 
marked deviation between theory and numerical results. That said we chose not to tune 


the velocity ratio appearing in equation (4.41 as (i) this is too speculative and (ii) such 
adjustment is certainly beyond the limits of our analogy. 

It is interesting to note that for £ smaller than /l, Mirel’s self-si milar variable r) m 
tends to Re 1 // = Z/Vt, so that the vertical structure for solution (4.51 now simply 
involves e rf ^ Re 1 //^. Actually, this solution is merely the purely diffusive solution of 
equation ( |4.2| ) for constant forcing ( u e constant in time). 

From the previous results we may extract several quantities, such as the displacement 
thickness or the locus of iso-velocities. The displacement thickness <5i can readily be 
estimated with f = erf as: 


p+oo 


5i = 


/Re . 


l_^dZ=-* 

u e J yWRe 


t — f 2 / 3 . 


(4.7) 


Similarly, isolines for the velocity can be extracted both for Qerris computations and for 
boundary layer theory. Figure [20] provides with a qualitative comparison between theory 
and numerical results, and it can be remarked that the overall prediction is more than 

just qualitative. _ 

Though the velocity ratio ( 2 r/n/\ / 3 t — f 2 )/(///2/t) (counterpart of Uoo/U s in equa¬ 


tion (4.41) is infinite near the shock, we note that the agreement between numerical and 


theoretical solutions is actually surprisingly good. We conjecture that whenever this ratio 
decreases to a value lower than one, i.e. near the centre and for large times where this 
ratio behaves as (4r)/(37r) so tends to 0, the error function approximation emerges as 
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Figure 19. Inner-boundary layer radial velocity profiles at different locations £ : 0.125, 0.25, 
0.5, 1 and 1.5. Blue solid lines correspond to numerical solutions obtained with Qtrris at 
t = 5 x 10 -3 ,10~ 2 ,5 x 10~ 2 and 10” 1 for Re = 5000 and We = 250 and represented in the 
self-similar space. Note that velocities are rescaled by their ma ximum value. The red dashed 
lines stand for the theoretical composite solution (equation ( |4.6| ) blending the self-similar vis¬ 
cous boundary layer solution with the self-similar Wagner inviscid solution fo r imp act. The 
composite solution is also rescaled by the edge velocity u e (f) given by equation (|3.30|). 


the solution of equation (4.41. Eventually we remark that an in-depth analysis of these 


phenomena demands a more involved description for the boundary layer (such as the 
Interactive Boundary Layer theory, see e.g. Lagree|20l0 ) and/or a deeper analysis of the 
Wagner region ( Qliver|[2002 Korobkin||2007 Qliver| 2007). 


4.3. Estimation of the shear stress and the total drag 
With this boundary layer solution, we are now in a position to provide with an estimation 
of the wall shear stress f = ^-dU r /dZ\^ , i.e. the viscous component of the stress which 


has been disregarded so far. And indeed this quantity is of paramount importance as far 
as raindrop-induced erosion of erodible beds is concerned ( Ellison|1945 Rein 1993; Lagree 
^003; Leguedois et al. 20051. Upon using equation (|4.5[) (with the same erf approximation 
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Figure 20. Left: Isolines for the radial velocity u r extracted from the numerical simulations. 
Right: Theoretical isolines for the radial velocity given by the composite expansion tt“ mp . Note 
that the transverse scale has here been stretched to visualize the boundary layer. 


for function f as before), we readily obtain: 


f(r, t) 


2 V3r 

7rlRe 1 ^ 2 (3t — f 2 ) 


(4.8) 


This theoretical prediction is confronted Fig. [21] with numerical profiles for the shear 
stress extracted from Qmis computations, and is shown to nicely agree with observations. 
From this local distribution for the stress we may infer the total drag induced with a 
drop impact, by integration over the wet area: 


n y/3t 

f(t) r dr d$. (4.9) 

Unfortunately this integral diverges due of the 1/x singularity developing in the near 
contact line region, and visible from Fig. [2l] left. Such singularities are usually a signature 
of an additional physics in the diverging region, not taken into account in the model. 
And indeed, Fig. [21] right reveals that the calculated shear stress significantly deviates 
from the theoretical prediction at some small distance A from the contact line position 
to reach a maximum value. Now integrating the local shear stress up to f = V3t — A, 
where A is this small cut-off length, we can provide an estimation for the drag at leading 
order in log(A): 


m = 3 VS (y 21 ° g (d) - 4 + log(12) )' (410) 

Upon noting that this quantity can be dimensionalised with pU 2 R 2 , the expression for 
the total drag in dimensioned variables follows: 


D(t) = —^p^U 2 RVt 
Vtt 



- 4 + log (12) 


) 


(4.11) 


Noticeably, the departure from the theoretical prediction pinpointed out in Fig. |21| right 
seems to occur at a precise location in self-similar variables, therefore suggesting a vl time 
dependence for A. From the numerical computations the value of A/Vt can be estimated 
to be around 0.03. Note that this is obviously a crude estimation, which nonetheless allows 
to propose the following estimate for the impact-induced drag: 


D{t) ~ 10 .7^p^U 2 RVt. 


(4.12) 


To further refine this prediction, the true nature of the cut-off length A needs to be 
clearly identified. Several candidates for governing this quantity naturally emerge, with 
for example the viscous 1/Re regularisation length in the vicinity of the contact line 
region or the inertial matching with the Wagner inner layer of typical size (d(t)/R) 2 ). 
This requires further investigation. 
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Figure 21. Left: Numerical and theoretical shear stress distribution underneath the drop, repre¬ 
sented in the self-similar space (the numerical data are taken at times i = 5xl0~ 4 , 10 -2 , 5xlCF 2 
and 10 _1 ). Right: Same data represented as a function of the distance from the contact line (log 
plot). This representation reveals a cut-off distance A from which the l/x singularity is screened. 
Importantly the numerical mesh size has been chosen to be small enough (Ax = 5 x 10 -4 ) to 
ensure the resolution of the fine-scale motion in the vicinity of the contact line. 
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5. Further comments and conclusion 

Capillary phenomena as well as possible aerodynamic effects from the surrounding gas 
have been disregarded so far. In this last part we shall estimate their influence on impact 
and discuss natural extensions of the present work. A summary and general conclusion 
then follow in §5.5| 


5.1. Influence of capillary phenomena 

The Weber number provides with a global measure of the ratio of available kinetic en¬ 
ergy to surface energy. For low values of this number (with respect to unity), drops 
typically bounce (Richard & Qucrc 20001 while preserving their shape or gently spread 
( Pasandideh-Fard et «/.|1996 1 according to the wetting properties of the underlying sub¬ 
strate. Note that even in this regime of droplet deposition, fast phenomena associated 
with the imbalance of surface stresses can set in ( Stebnovskii| [T979). When the initial ki¬ 
netic energy of the drop can no more be neglected in front of the surface energy (We ~ 1), 
a surface wavefield starts to develop on the drop, shaping it into a characteristic liquid 
pyramid or torus (Renardy et al. 2003). For even higher values of the Weber number, such 
as the inertial limit We 1 investigated in the present paper, we do not expect capillary 
phenomena to have a significant influence on a global scale, but locally surface tension 
can still play a dominant role. For example, the high-curvature turnaround region at the 
lamella root is typically a place where capillarity presumably plays an important role. 
But due to scale separation, this region is invisible at our level of description. Indeed, 
in classic impact analyses see e.g. Oliver (20021, the typical extent of this intermediate 
Wagner region associated with highly curved interfaces is found to be 0(e 2 ), to be com¬ 
pared both with the O(e) size of the main impact region considered throughout this paper 
(see I 2.31 and with the 0(e 3 ) thickness of the lamella. In the framework of our first-order 
theory, we therefore do not anticipate appreciable deviations stemming from this zone. 
Conversely, for a correct description of the ejected liquid sheet feeding conditions and of 
the pressure fall-off near the lamella root reported Fig. |15[ an accurate representation of 
this matching region appears mandatory. 

The contact line is an other region where marked effects from capillarity are to be ex¬ 
pected. Drop impact is characterized with fast motions near the contact line. This violent 
dynamic wetting phenomenon can arguably bring about issues in our description of im¬ 
pact. Actually Blake et al. (1999) demonstrated that nonlocal hydrodynamics could play 
a significant role in the dynamic contact angle selection. Based on experimental data, 
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|Blake et al.\ further put forward the possible ‘mutual interdependence’ between the phe¬ 
nomena in the near contact line region and the far-held hydrodyn amics. This complex 
interplay was further confirmed in the context of drop impact by Sikalo et al. (2005), 


but especially for the late receding phase. Interestingly, these authors demonstrated that 
the early evolution of the contact angle was quite insensitive to the experimental con¬ 
ditions and fairly well captured by the contact angle of a truncated sphere. This nice 
agreement certainly advocates for a predominance of inertial effects over capillary cor¬ 
rections emanating from the dynamic contact line, at least in the early stage of impact. 
And indeed, remembering that shortly after impact the fluid motion in the contact area 
is essentially vertical, it appears likely that the point of contact can be determined with 
mere inertial arguments. In our simulations, dynamic effects have been disregarded in 
the description of the contact angle, which has been set to the constant value 7t/2. The 
agreement between our simulations and the purely inertial theory is again an indication 
of the unimportance of dynamic wetting. It might further be interesting to note that the 
surface energy gained by wetting the solid is of the order of 1/We when rescaled by the 
initial kinetic energy. Again this heuristically rules out any leading effect from wetting 
in the short-term dynamics. This ratio evolves with time though, and ultimately wet¬ 
ting phenomena become dominant, as evidenced by the late f 1 / 10 spreading behaviour 


consistent with Tanner’s law in the experiments of Rioboo et al. (2002). 


5.2. Influence of ambient air 

For about a decade or so, there has been an increasing realization of the role played by 
surrounding air in liquid impact in general, and drop impact in particular. Following key 


experiments performed by Xu et al. (2005) on air-induced splash triggering, a number of 


studies have focused on the events preluding liquid sheet ejection. The first significant 
effect of surrounding air is to impart a dimple-like deformation in the bottommost region 


of the drop (see experimental observations of the dimple obtained by Thoroddsen et al. 


2005 and X-ray ultra-fast imagery of the complex dynamics of this air pocket by Lee 


et al. 2012). Smith et al. (2003) first depicted theoretically this process by coupling 


lubrication in the squeezed air film and potential flow inside the drop. These authors 
notably evidenced the presence of off-axis pressure peaks. While more recent studies 
raised doubt about the link between this dimple formation and splash triggering per se - 
that might merely be a secondary independent consequence of the presence of surrounding 
gas ( Duchemin fe Josserand|2011 1, this gas pocket is nonetheless formed over timescales 
and lengthscales overlapping that of the phenomenon reported in the present paper (Mani 
et al. 2 0101 ). It is therefore legitimate to question about the impact of this air entrapment 


phenomenon on our results. 

In order to investigate these effects, we performed a simulation of a liquid drop approach¬ 
ing a solid substrate, deforming as a result of lubrication pressure rise in the film, and 
finally impacting the substrate. Figure [22] represents the time evolution of the pressure 
exerted on the support at the axis. Considering that air delays the moment of impact 
( Mani et aZ.|2010 ), we introduce a time shift fimpact corresponding to the moment where 
the drop and the solid are only a grid cell apart. Interestingly our results reveal that 
the pressure at the origin (measuring now the entrapped bubble pressure) is fairly well 
captured by relation (3.371 after replacing t with the true time from impact t — fi mpa ct, 
that is: 


/q 

p(0, 0,t) = - {t /impact) 


(5.1) 


This agreement between our prediction and a simulation incorporating air entrapment 
effects not only validates and extends our results beyond the initial scope of Wagner 
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Figure 22. Left: Time evolution of the pressure as measured under an impacting drop with 
air-induced dimple formation (i.e. bubble entrapment) taken into account. The red trace moni¬ 
tors the pressure at the origin. The physical parameters of this Qerris simulation are Re = 5000 
and We = 250. The superimposed black dashed line corresponds to the theoretical solution 
p(0,0, t) = ^(t — (impact) - 5 delayed from (i mpa ct, where (impact is the real impact time (in the 
numerical simulations, (impact corresponds to the time at which liquid and solid are just one grid 
cell away). Right: Close-up of the bottommost point of the drop in the numerical simulation (at 
i = 1.54 x 10“ ). The position of the interface is materialised with a blue line. The colormap 
illustrates the distribution of the pressure field within the drop and in the gas layer. Noteworthy 
enough the isopressure lines seamlessly cross the interface, revealing the transparency of the 
dimple to pressure. Note that for the sake of clarity, the vertical scale has here been magnified 
by a factor 22. 


impact theory (disregarding air effects), but also suggests that the results of the present 
manuscript correspond to the far-held behaviour of an impacting drop in presence of 
surrounding gas. This observation outlines the appealing prospect of describing both the 
dimple geometry and associated dynamical fields by analytical means. 

5.3. Main results 

In this paper, the short-term dynamics of a drop impacting a rigid substrate has been 
elucidated. A self-similar solution for the impact-induced how has in particular been 
unraveled and matched to a self-similar viscous boundary layer. This solution has been 
intensively validated with numerical Qerris computations, and this constant cross-testing 
between asymptotic theory and multiphase adaptive how simulations is one of the key 
feature of the present approach. In the course of this investigation, several important 
results have been substantiated. These results allow both for a simple yet accurate qual¬ 
itative depiction of drop impact along with an in-depth quantitative understanding of 
this phenomenon. These key results are summarised in the following: 

• A fundamental analogy between the water entry of a solid object (Wagner’s original 
problem) and drop impact exists, 

• During the earliest moments post-impact, the contact line follows a tank-treading 
motion. There is in particular no contact line sweeping motion, 

• The impact-induced how is concentrated in the contact zone, and the far-held merely 
corresponds to an undisturbed rigid-body motion reducing to a global free-hight at ve¬ 
locity U. There is no global or large-scale drop deformation during impact, 

• The position of the contact line is given by the simple relation d(t) = V3RUt. 
Though simple, this locus does not correspond to the cut radius of a truncated sphere, 

• The wet footprint extent of the drop dictates the size of the impact-induced per¬ 
turbed how, 
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(r = 0, z == 0(d(t)) > 0 + ) 


(0 < r < d(t) = \/3RUt, z = 0 + ) 




u r (0, z,t ) = 0 


u~(0,z, t) = ™ arctan ^ V3URt ^ 


— arctan 


V3U 3 Rt ~ _ TT 
3URt+z 2 u 


u z (r, z, t) = 0 


( d(t ) » r, d(t) » z > 0 + ) 


(0 < r < d(t) = VSRUt, z = 0(Re~ 1/2 d(t)) 




ti(URt—ZZ-) ^ 


P(0,0,t) = £^ 


3 1 

2fl2 



F(t) = 6V3puiRiVt D(t) ~ 10.7pip? U 2 RVt 


Table 2. Summary of the main results of the paper in dimensioned form. The top part of the 
table refers to ideal fluid results (left: closed-form results along the axis of symmetry, right: 
along the substrate). The left middle part sums up inviscid stagnation point results, and the 
right middle part summarises the viscous boundary layer results. Observables such as the net 
normal force F(t) and tangential force D(t) are reminded as well in the bottom part of the table. 


• There is a consistent analogy between the impact-induced flow within the drop and 
the flow induced by a flat rising expanding disk (Lamb’s analogy), 

• The impact pressure is to be associated with the unsteady Bernoulli contribution 
—dt4>- It cannot be inferred from usual inertial steady contribution — pU 2 , 

• As a corollary to the previous point, the impact pressure is extremal at the contact 
line. It is not maximal at the stagnation point, 

• A full three-dimensional self-similar solution for the impact-induced flow of an in¬ 
viscid drop exists and matches quantitatively realistic numerical data on drop impact, 

• Analytical solutions for this flow have been presented in integral forms (with some 
explicit closed-form expressions along some particular locations), see table [573] (in dimen¬ 
sional form), 

• An original inviscid stagnation point structure with an unexpected r/ \ft. slip velocity 
develops in the vicinity of the origin. The velocity field structure markedly differs from 
the classic r/t prediction occurring for later times, 

• An approximate self-similar solution for the viscous boundary layer seamlessly match 
with the inviscid impact-flow (analogy with Mirels shockwave problem), 

• Self-similar variables have the same structure z/\/t both in the outer region and in 
the boundary layer, 

• From the knowledge of the distribution of the dynamical fields across the wet area, 
the expressions for the normal and tangential total force on the substrate are provided, 

• The asymptotic solution is found to be numerically valid over several decades in 
time. This solution was found to be insensitive to air-induced dimple formation, 

• For times of order one, the present results remain at least qualitative. 
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The results obtained in the present manuscript offer several appealing prospects. On the 
role of surrounding air first, the last results of ( }5.2| support the idea that the dimple 
geometry and characteristics could be derived analytically, with a far-field corresponding 
to the here presented flow. The question of the role of air in suppressing the pressure 
divergence (Josserand et al. 2015), or in altering the shear stress distribution is also a 
central point for a correct description of a single drop impact action. Similarly, the ques¬ 
tion of the scope of the present results for impact on a liquid film or on a soft/erodible 
substrate is also of interest. Indeed, the analytical and numerical toolboxes developed 
here might well be transposed to other rheologies (such as Bingham, see |Staron et al. 


(20131 or granular media, ( Lagree et q?.||20lT )). While the inverse-square-root singular¬ 
ity of the pressure is integrable and yields a finite normal force, the sharper singularity 
observed in the viscous shear stress distribution leads to a divergence for the total drag. 
From the simulations, it appears that the singularity is screened over a lengthscale A, 
but the precise nature of this cut-off scale is still unclear. As discussed in both vis¬ 
cous effects and inertial ones are candidates to explain this regularisation. While the 
rudimentary description for the viscous boundary layer proposed in the manuscript cer¬ 
tainly necessitates a refined analysis, the question of inertial effects at the Wagner region 
scale is also worth studying. Not only an in-depth investigation of such effects might 
provide an explanation for the regularisation of the shear stress, but it shall shed light 
over the onset of lamella formation, which still conceals mysteries. Finally the process 
responsible for the loss of self-similarity observed at intermediate times is still uncertain: 
confinement effects arising when the impact-flow lengthscale overlaps with the drop ex¬ 
tent, eventual deceleration in the far-field region or contact line geometrical departures 
from the square-root law are all a priori legitimate to explain the final pressure fall-off, 
and certainly needs further investigation. 


5.5. To conclude 

Within the numerous limits carefully drawn along this paper, a consistent asymptotic 
description of the dynamics and geometry of drop impacting a solid surface has been 
proposed. The results may simply be summarised through three analogies: Wagner water 
entry (drop impact being the dual of this problem), Lamb’s disk winding flow (that 
accurately represents the flow induced with the impact) and Mirels shockwave-induced 
boundary layer (remarkably capturing the boundary layer developing in the contact 
line’s trail). The original strategy developed throughout this paper has been to validate 
those three analogies through a constant confrontation between numerical simulations 
and asymptotic analysis. Our study revealed that very powerful state-of-the-art adaptive 
codes now allow to probe all the dynamic features of realistic violent events such as 
drop impact, but in the meantime, it also emphasised again how powerful and useful 
asymptotic analysis is in providing an in-depth understanding of such phenomena and in 
uncloaking the raw data delivered by the code. Finally our study brought to light some 
interesting features and observables (such as the particular stagnation point structure, 
pressure distribution, contact line motion, viscous total drag force) never observed to date 
neither in simulations nor in experiments. This certainly arouses the exciting prospect 
of their unveiling in future experimental studies. 
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Figure 23. Typical mesh structure refined adaptively by Qerris flow solver during a simulation. 


6. Appendix 

6.1. Qerris flow solver 

All the numerical simulations were performed with the open-source code Qerris (freely 
downloadable at http://gfs.sourceforge.net - see also Popinct 2003, 2009; Lagree 


et al. 2011 for details). Qerris is a solver of the incompressible Navier-Stokes equations 


taking into account multiple phases and surface tension. The code makes use of a finite- 
volume approach and of a Volume-of-Fluid (VoF) method for an accurate description 
of the transport of the interfaces between two-phase flows. It also features an adaptive 
mesh refinement procedure allowing for both a precise description of flows with large 
scale separation and a reduction of computational costs. Typically in our simulations the 
finest grid is chosen to be concentrated along free surfaces and within the contact zone to 
fully capture the features of the pressure field and of the boundary layers (see Fig. [23). In 
these areas the corresponding local resolution usually corresponds to 4096x4096 but can 
reach local density as high as 32768x32768 if needed (examples being Fig. [l] or Fig. [2). 

The simulations carried out in this study all correspond to the impact of a water drop 
in air with a Reynolds number of 5000 and a Weber number of 250. The computations 
were performed in an axisymmetric configuration. We emphasize that both liquid and 
air motions were computed with Qerris, but, to be consistent with the post-impact theory 
developed in this paper the simulations disregarded air cushioning and dimple formation 
(except explicitly specified, see (5.2). To avoid dimple formation in this multiphase flow 
simulation, the initial configuration is set to a slightly truncated liquid sphere already 
touching the solid surface. The liquid is initialised with a constant downward velocity. 
The initial sphere penetration ro = 10 -4 is at most one grid cell deep (for example the 
grid spacing is Ax ~ 5 x 10 -4 for 4096x4096 simulations). Finally a no-slip boundary 
condition is enforced at the substrate level. The reliability of the results has been thor¬ 
oughly checked with a convergence study on the refinement level, and with particular 
attention paid for the pressure field and the position of the contact line convergence. 
Fig. [24] proposes a comparison between the evolution of the pressure field measured 
at the origin for two levels of resolution (2048x2048 and 4096x4096). For both cases 
the numerical solut ion quickly converges to the theoretical solution j?(0,0,t) = ^ 

(see equation (3.37)) around t = 5 x 10 -3 and leaves the self-similar regime at around 
t = 6 x 10“ 1 . We remark that after a transient period, both numerical solutions give con¬ 
sistent information and collapse onto the theoretical solution over almost three decades. 
Note that the occurrence of sporadic glitches in the numerical solution (see e.g. Figs 
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Figure 24. Comparison between the pressure measured in the Qerris simulations at the origin 
for two different maximum level of mesh refinement (red dots) and the theoretical prediction 
(black dashed line). The left panel corresponds to a simulation where the maximal grid density 
is 2048x2048, and the right panel to a simulation where the maximal density is 4096x4096 (in 
both cases the physical parameters are Re = 5000 and We = 250). The analytical solution of 


the pressure is given by p(0, 0, f) 


_ V3 


t 2 , see equation (3.37l. After a short transient, both 


simulations quickly reach the same self-similar asymptotic regime. 



t 


Figure 25. Close-up of Fig. |24| right. Note that the numerical evolution for the pressure is 
slightly above (about 7 %) the analytical solution. 


17 24 or 25) are to be related with the classic difficulty of computing the pressure in 


projection methods, such as the one implemented in Qerris (Brown et al. 2001 Popinet 


2003). We finally remark that an error of ca. 7 % between the numerical prediction for 


the pressure at the origin and the theoretical prediction was consistently noted in our 
simulations (see Fig. 25). The nature of this discrepancy is uncertain though, and might 
either be related to the aforementionned numerical difficulties in computing the pressure 
or to the limits of our first-order asymptotic description for drop impact. 


6.2. gerris parameter file 

A minimal Qerris parameter file allowing to reproduce the results presented in this manuscript 
is provided here for the reader’s convenience (Note for referees and Editors: this section 
will likely be suppressed in the final version for this script will be on a dedicated example 
page in Qerris website): 

Define tO le-4 
Define Re 5000 
Define ReAir 277778 
Define We 250 

Define VAR(T,min,max) (min + CLAMP(T,0,1)*(max - min)) 

Define RH0_EAU 1000. 

Define RH0_AIR 1. 

Define RH0(T) VAR(T, RH0_AIR/RH0_EAU, RH0_EAU/RH0_EAU) 
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4 3 GfsAxi GfsBox GfsGEdge {y=0.5 x=0.5} { 

Time { t=tO end = 2.5 dtmax = le-4 } 

PhysicalParams { L = 2 } 

VariableTracerVOFHeight T 
VariableFiltered T1 T 1 
VariableCurvature K T Kmax 
SourceTension T (l./We) K 
InitFraction T ({ 

double goutte = -(x+tO-1.)*(x+tO-1.) - (y)*(y) + 1.*1.; 

return (goutte); 

}) 

Init {}{ U=-1*(T) V=0} 

# Initial refining 

RefineSurface 12 ( -(x+tO-1.)*(x+tO-1.) - (y)*(y) + l.*l. ) 

Refine ( ((x>0.)&&(x<0.3)&&(y>0.)&&(y<4.)) ? 12 : 4 ) 

# Refining on tracer gradient 

AdaptGradient { istep = 1 } { maxlevel = 12 cmax = le-2 } T 

# Constant value of mesh’s level inside the bulk 

AdaptFunction { istep = 1 } { 

maxlevel = ( ((x>0.)&&(x<0.3)&&(y>0.)&&(y<4.)) ? 12 : 4 ) 
cmax = le-2 } (T==l) 

RemoveDroplets { istep = 1 } T -1 
PhysicalParams { alpha = 1./RH0(T1) } 

SourceViscosity { } (T*(l./Re) + (1. - T)*(1./ReAir)) {beta = 1} 

OutputTime { istep = 1 } stderr 

OutputTiming { step = 0.0001 } stderr 

Outputsimulation { istep = 1 } stdout 

OutputProjectionStats { istep = 1 } stderr 

} 

#1 

GfsBox { 

bottom =Boundary{} 
left = Boundary{ 

BcDirichlet U 0 
BcDirichlet VO} 

} 

#2 

GfsBox { 

right = Boundary{ 

BcDirichlet P 0 
BcNeumann U 0 
BcNeumann VO} 
bottom = Boundary{ } 
top = Boundary{ 

BcDirichlet P 0 
BcNeumann U 0 
BcNeumann VO} 

} 

#3 

GfsBox { 

left = Boundary{ 

BcDirichlet U 0 
BcDirichlet VO} 
right = Boundary{ 

BcDirichlet P 0 
BcNeumann U 0 
BcNeumann VO} 


#4 


} 
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GfsBox { 

right = Boundary! 

BcDirichlet P 0 
BcNeumann U 0 

BcNeumann VO} 
top = Boundary! 

BcDirichlet P 0 
BcNeumann U 0 
BcNeumann VO } 
left = Boundary! 

BcDirichlet U 0 
BcDirichlet VO} 

} 

1 2 right 
1 3 top 
3 4 top 
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